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ABSTRACT 

Using density functional molecular dynamics simulations, we determine the equation of state 
for hydrogen-helium mixtures spanning density-temperature conditions typical of giant planet 
interiors, ^ 0.2— 9gcm~^ and 1000— 80 000 K for a typical helium mass fraction of 0.245. In 
addition to computing internal energy and pressure, we determine the entropy using an ab initio 
thermodynamic integration technique. A comprehensive equation of state (EOS) table with 391 
density-temperature points is constructed and the results are presented in form of two-dimensional 
free energy fit for interpolation. Deviations between our ab initio EOS and the semi-analytical 
EOS model by Saumon and Chabrier are analyzed in detail, and we use the results for initial 
revision of the inferred thermal state of giant planets with known values for mass and radius. 
Changes are most pronounced for planets in the Jupiter mass range and below. We present a 
revision to the mass-radius relationship which makes the hottest exoplanets increase in radius by 
0.2 Jupiter radii at fixed entropy and for masses greater than ~ 0.5 Jupiter mass. This change 
is large enough to have possible implications for some discrepant "inflated giant exoplanets" . 

Subject headings: equation of state, hydrogen-helium mixtures, ab initio simulations, giant planets, 
extrasolar planets 



Introduction 



Th e semi-analytical model bv lSaumon fc Chabrier 
(|l992l ) (SC) for the EOS of hydroge n and its exten 



sion to hydrogen-helium mixtures (jSaumon et al 



19951) were very successful and have been used 
in numerous calculations for the interiors of gi- 
ant planets. However, with the development of 
ab initio computer simulation techniques many 
uncontrolled approximations can now be avoided, 
simplifications inherent to analytical EOS models 
and severely limiting their predictive capabilities 
in the regime of high density and low temperature 
where interactions between particles are strong. 
Relying solely on analytical methods, it is difficult 
to determine the ionization state of the differ- 



ent chemical species that are present in the dense 
fluid. 

Ab initio simulations allow one to study a fully 
interacting system of particles and to determine its 
properties by deriving the electronic states explic- 
itly for every configuration of nuclei. No parame- 
ters are adjusted to match experimental data, but 
ab initio simulations still rely on approximations 
to solve the Schrodinger equation. However, they 
are not specific to the particular material nor the 
pressure-temperature conditions under considera- 
tion. 

In this paper, we rely on density functional 
molecular dynamics (DFT-MD) simulations that 
have been employed before to study hydro- 
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Collins et al. 



Mor ales et al 
2012t" 



20inl: 



gen (iLenoskv et alj 120001: iMilitzer et all 12001 
Desiarlais'20 03tlBoney et al.|200 4':'Nettelmann et_^al 

2008: 



Caillabct e t a l. 20 IJ 
Nettelmann et alJ |2012|). he- 



proxi mation ( Saumon et al.ll995t iNettelmann et al 



2008[ ). we computed the EOS over a wide range 



lium (iMilitzeij 120061: IStixrude fc Jeanlozj |2008l; 

Militzer 12009 ) and hydrogen - helium mixtures (.Vorberger 



2011I ). While the computation of the pressure 
and the internal energy is straightforward from 
DFT-MD simulations, the entropy is not directly 
accessible. However, an accurate knowledge of 
the entropy of hydrogen-helium mixtures at high 
pressure is of crucial importance for the deter- 
mination of the temperature profile, the density, 
and the thermal energy budget in the interior of a 
giant planet. In 2008, two groups constructed 
Jupiter int erior mode l s from DFT-MD simu- 
lations (iMilitzer et al.l l2008t INettelmann et al 



20081 ). While the derived pressures and inter- 
nal energies can be considered to be more reliable 
than those predicted by the SC model, both pa- 
pers predicted very d ifferent interior temperatur e 
profiles for Jupiter (jMilitzer fc HubbardI [iooi). 
Using ab initio thermodynamic i ntegration tech - 



niques (TDI), we recently showed (lMilitzeill2013f) . 
that the work bv INettelmann et al.l (|2008[) " overes- 
timated the temperature at Jupiter's core-mantle 
boundary (CMB) by 3050 K (19%) while we un- 
derest imated it by 2870 K (18%) in IMilitzer et al.l 
( 20081 ). The revised temperature for the Jupiter's 
CMB is 16 150 K and the corrections to the SC 
EOS model are in fact only — 350K. 

At conditions of Jupiter's CMB, hydrogen is 
metallic and characterized by a high degree of elec- 
tronic degeneracy. Such a degenerate state is de- 
scribed rather well by the the SC model. How- 
ever, when we applied the TDI technique to ex- 
plicitly determine the entropy over a wide range 
of pressure-temperature conditions, we identified 
a number of discrepancies between the DFT-MD 
results and the SC model. Near the molecular-to- 
metallic transition, our simulations predict a sig- 
nificant shift of the adiabat towards higher densi- 
ties. At high temperature, where electronic exci- 
tations matter, our computed entropies are higher 
than those of the SC model. We also do not per- 
fectly reproduce the SC entropies in the molecular 
regime at low density. 

Rather than providing a separate hydrogen and 
helium EOS and relying on the linear mixing ap- 



of density-temperature conditions for a represen- 
tative mixing ratio of A^ho =18 helium atoms 
in A h = 220 hydrogen atoms, corresponding 
ettffll.a helium mass fraction of y =0.245, which 
is close to the solar value. This means that 
the nonideal mixing effect s are f ully incorpo- 



rated. In Vorberg er et al.l (|2007bl ). we showed 
for example that the presence of helium makes 
the hydrogen molecules more stable and re- 
duces the dissociation fraction at given pres- 
sure and temperature. Even if other mixing ra- 
tios b ecome of interest, as the r e sult of helium 
rain ( Stevenson &: Salpeted 1977| : Morales et al 



2009; Wilson fc Militzer 



2010t 



McMahon et al 



|2012[ ). one is still better off by starting from an 
EOS for a typical hydrogen- helium mixture and 
then perturbing the mixing ratio by a compara- 
tively small amount. Increasing or decreasing the 
helium fraction requires knowledge of a helium 
or hydrogen EOS, respectively. For the helium 
EOS, we recommend our first-principles compu- 
tation (|Militzeij [ioogh because it provides simu- 
lation data points for the pressure, P, internal 
energy, E, Helmholtz free energy, F, and entropy, 
S, over a wide parameter range and a thermody- 
namically consistent free energy fit for interpola- 
tion. For available hydro gen EOS work, we refer 
to the recent review by iMcMahon et al. ( 20121 ) 



but there has also been a considerable theoreti- 
cal effort compute th e hydrogen EOS with semi- 
analytical techniques (iDharma-wardana fc Perrot 

12002: K raeft et al . 2002': 'Rogers fc Navfonov'2002t 

Sa fa fc Pfennigeril200& : .Ebehng et al.il2012c lAlastuev fc Ballenegger 



20121 ). If the perturbation in the helium fraction 
is sufficiently small, one may use the SC EOS for 
simplicity. 

2. Ab Initio Simulations 

We base our a b initio entrop y calculations on 
our recent article ( Militzeiil2013l ) where we showed 



how the TDI technique can be extended to study 
molecular hydrogen and how it can be applied effi- 
ciently to determine the entropy at high tempera- 
ture where electronic excitations matter. The TDI 
technique allows one to determine the difference in 
the Helmholtz free energy between two interact- 
ing many- body systems at fixed density and tem- 
perature ( Morales et al. 20091 : Wilson fc Militzer 



2 



2010l l2012allbl : iMcMahon et al.ll2012h . We apply 
this method to determine the free energy differ- 
ence between the DFT simulations and a system 
of classical forces that we construct: 



^DFT 



dX {Vks - 



(1) 



The angle brackets represent an average over tra- 
jectories governed by forces that are derived from 
a hybrid potential energy function, Vx = XVks + 
(1 — A)yci- Vci is the potential energy of the 
classi cal system and Vks is the Kohn-Sham en- 
ergy (iKohn fc ShamI [19651 . The presence of elec- 
tronic excitations leads to an intrinsic contribu- 
tion to the entropy and aff ects the forces on the 
nuclei (|de Wiis et al.l 119981 ) that need to he de - 
rived from the Mermin free energy ()Merminlll9"65l) . 
^ = ^KS — We combined both contributions 

into t he following ex pression for the ab initio en- 
tropy (|Militzeij|2013l ): 



TS = {Vks) + (/^io„) ~ f dX {n - V,i)^ - Fd. 

JQ 

^ (2) 

(^s) includes contributions from partially occu- 
pied excited states. The A integration was per- 
formed using five independent MD simulations 
with A equally spaced between and 1. To make 
this integration process efficient, we construct the 
pair potentials of the classical syste m to match the 



DFT forces as closely as possible (jlzvekov et al 



2003I ). The computation of classical free energy is 
performed with Monte Carlo methods by thermo- 
dynamic integration to an system of noninteract- 
ing particles. 

All simulations were performed w ith the VASP 
code ( Kresse fc Furthmiiller 1996[ ) with pseu- 
dopotentials of th e projector-augmented wave 
type (jBlochll 119941) and a plane wave basis set 
cutoff of at least 1000 eV. The Perdew-Burke- 



Ernz erhof exchange-correlation functional (jPerdew et al 
19961 ) was used throughout, but it has been shown 
recently that simulations based on the local den- 
sity approximation yielded very similar r esults 
for Jupiter's deep interior (|Militzeii 120131 ). In 
the same article, we also performed a combined 
finite-size and k point analysis that demonstrated 
that simulations with 256 electrons and the zone- 
average point k — are sufficiently accu- 
rate. All results that we report in this article were 



thus obtained with 220 hydrogen and 18 helium 
atoms in periodic boundary conditions. 

We used a MD time step 0.2 fs, except for tem- 
perature of 50 000 K and above where we used 
a time step of 0.1 fs to accurately capture the 
more rapid collisions between particles at elevated 
temperatures. All standard DFT-MD simulations 
that we performed to determine P and E were 
2.0 ps long, except at the highest temperatures, 
where 1.0 ps were found to be sufficient because 
the auto-correlation times are short and the er- 
ror bars are small. All simulations were initialized 
with positions and velocity vectors from converged 
MD simulations at nearby densities and tempera- 
tures. This allowed us to run the TDI simulations 
for only 0.5 ps at each A point. 

We also adjusted the number of orbitals in the 
calculations to accommodate the partial occupa- 
tion of excited e lectronic states according to Mer- 
min functional ( Mermin 19651 ). The number of 



orbitals was increased until the error in the inte- 
gral of the Fermi function was reduced to less than 
10~^. This required many orbitals at high tem- 
perature and low density. Up to 816 were used, a 
significant increase in the computational cost over 
the 128 needed for ground state calculations. This 
is the primary reason why we omitted simulations 
that would lead to entropy values above approxi- 
mately 12.5 fcf,/el. The regime of higher temper- 
atures can be studied much more efficiently with 
path integral Monte Carlo (PIMC) simulations 
because the computational cost of this alterna- 
tive first-principles simulation technique scales like 



droeen (IPierleoni et al.' 'l994'; 


Magro et al. 


1996 


Militzer et al. 1999; ,Militzer & Ceoerlev 


200C 


Militzer & GrahamI 20061 Hu et al. 12010, l201ll) 



hefium CMilitzer 2006, 2009^, and hydrogen- 
helium mixtures (Militzer. 2005,) at high pressure 
and temperature and most rece ntly also to study 
the EOS of carbon and water (Driver fc Militzer 



20121 ) 



3. Equation of State Results 

We report the computed equation of state in 
the form of a table, a series of figures, and in ana- 
lytical form as two-dimensional fit of the free en- 
ergy. In table [TJ we provide the thermodynamic 
functions that directly follow from analysis of the 
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Fig. 1. — Temperature-density conditions of DFT- 
MD simulations. The circles indicate parameters 
where entropy and free energy have been calcu- 
lated in addition to the pressure and internal en- 
ergy. The lines show adiabats. The labels specify 
their entropy values in units of kj, per electron. 
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Fig. 2. — Internal energy per electron as function 
of temperature for four different densities given in 
terms of r^. Results from DFT-MD simulations 
are compared with the analytical SC model. 
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Fig. 3. — Internal energy per electron as function 
of density for seven different temperatures. Pre- 
dictions from DFT-MD simulations and from the 
analytical SC EOS model are compared. 
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DFT-MD trajectories. The pressure and internal 
energy were computed for 391 different density- 
temperature points (see Fig.[T]). Tlie 1-a error bars 
correspond to statistical uncertainty that arises 
from the finite length of the MD simulations. For 
131 points in table[Tl the thermodynamic integra- 
tion was performed with five A points and the free 
energy and entropy are reported in addition. Only 
counting the production runs that led to results 
in table 1, the total CPU time consumed for this 
project amounted to 850 000 core-hours on Intel 
Nehalem processors. This is equivalent to using 
100 cores for an entire year, which is a consider- 
able amount of computer time by today's stan- 
dards but will certainly become available to ev- 
eryone in the near future as computers with more 
and more cores are assembled. 

In figuresHHHElllEllHl andM we plot the in- 
ternal energy, pressure, Helmlioltz free energy, and 
entropy respectively as a function of temperature 
and density. Every circle corresponds to a particu- 
lar DFT-MD simulation listed in table [U without 
any interpolation being performed. The dashed 
lines are the results of the most common version 
of the analytical SC EOS model where the differ- 
ent thermodynamic functions have been smoothly 
interpolated across the molecular-to-metallic tran- 
sition in hydrogen. 

To accommodate the wide parameter range of 
our simulations, we plot the different thermody- 
namic functions on logarithmic scale. Since these 
functions depend strongly on density and temper- 
ature, we added a second panel where we removed 
most of this dependence by introducing a scale fac- 
tor equal to or T raised to some power. Here 
is the Wigner-Seitz radius that specifies the den- 
sity of system according to = V/N^ = n^^, 
while n is the number of electrons, TVe = {Nn + 
2A^Ho), per unit volume V. The mass density is 
given hy p = n(7VHmH + Nii^mac) / {Nu + 2Niic) 
where mu and tohc are the masses of the hydrogen 
and helium atoms. 

The rescaling of the ordinate makes it eas- 
ier to identify the deviations from the SC model 
while our simulation results can still be repro- 
duced easily. The ordinates are plotted in atomic 
units. Lengths including rs are given in Bohr radii 
(ao=5.29177209xl0~"m), energies in Hartrees 
(4.35974380x10-18 J) per electron (el.), and en- 
tropies are specified in units of kb per electron. 



where kb is Boltzmann's constant. 

Figure [2] shows a comparison between the inter- 
nal energies from DFT-MD simulations with the 
predictions of the SC EOS model. For a low den- 
sity of = 2.4, excellent agreement is found for a 
temperature range from 1000 to 20 000 K. Hydro- 
gen gradually changes from a molecular state to 
an atomic state in this temperature interval and, 
from the good agreement, one may conclude that 
the thermally activated dissociation of molecules 
is well described in the SC model. However, 
20 000 K, the SC model predicts an strong and ar- 
tificial increase in the internal energy that is the 
result of an inaccurate description of electronic 
exc itations. This devi a tion w as first identified 
by mnt zer fc CeperlevI (|200lh when predictions 
from the SC model were compared with PIMC 
simulations. Figure |3] shows that this deviation is 
present at 20 000 K for whole density range under 
consideration and extends to much higher temper- 
atures also. 

Figure [2] shows that the favorable agreement 
between DFT-MD results and SC predictions be- 
low 20 000 K continues to hold up to a density of 
rg = 1.6. When the internal energy is compared 
for a higher density of = 1.0 or 1.2 where hy- 
drogen is metallic, one finds that DFT-MD results 
and SC predictions are offset by a nearly constant 
amount. 

The internal energy curves of = 1.6 and 2.4 
appear to cross over in Fig. [2] at a temperature of 
27 000 K, which is consistently predicted by DFT- 
MD results and the SC model. Figure [3] shows 
that this is simply a consequence of internal en- 
ergy exhibiting a minimum when plotted at con- 
stant temperature as function of density. At high 
density, the internal energy sharply rises because 
of Pauli exclusion effects between the electrons. In 
the low density limit, the internal energy rises also 
because the ionization fraction increases as a result 
of the increased gain in entropy that is associated 
with electrons becoming free particles. 

In Fig. m we compare the pressure predicted 
from DFT-MD simulation with the SC model. At 
a high density of = 1.0 where the hydrogen- 
helium mixture is metallic, we find fairly good 
agreement over the entire temperature range. This 
implies that the deviation that we identified for 
the internal energy in this regime, varies slowly 
with density and does not significantly affect the 
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Fig. 4. — Pressure isochores computed with DFT- 
MD simulations are compared with the analytical 
SC model. 
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Fig. 5. — Pressure isotherms computed with DFT- 
MD simulations are compared with the analytical 
SC model. 
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pressure in the SC model. 

At a low density of Ts = 2.2 and 2.4, we found 
good agreement up to a temperature of 5000 K. 
At this temperature, we see a small decrease in 
slope in the DFT-MD data that is missing in the 
predictions of the SC model. We attribute this 
slope change to the dissociation of molecules in the 
DFT-MD simulations. At 20 000 K, the SC model 
predicts a significant decrease in slope which is not 
present in the DFT-MD data. This slope change 
in the SC predictions can again be attributed to an 
inaccurate description of ionization, which leads to 
deviations over the whole density range under con- 
sideration (Fig. [5|) . At an intermediate density of 
Ts = 1.6 close to the molecular-to-metallic transi- 
tion, we find that the SC model overestimates the 
pressure up to about 20 000 K and underestimates 
for higher temperatures. The deviations around 
100 GPa, 5000K, and r, = 1.6 (0.75 gcm^^) are 
of particular significance. The DFT-MD simula- 
tions predict pressures that are much lower than 
those of the SC model. This leads a significant 
departure in the resulting adiabats. Its implica- 
tion for the interiors of giant planets will later be 
further analyzed. 

In figures [6] and [3 the Helmholtz free energy 
from DFT-MD simulations and the SC model are 
compared. In general the agreement appears to be 
much better than for other thermodynamic func- 
tions that are derivatives of it. Still, one finds that 
the SC model overestimates the free energy in the 
metallic regime, mirroring the deviations that we 
have discussed for the internal energy. 

In figures [51 the entropies at different densi- 
ties are compared as a function of temperature. 
At a very high density of = 1.0, very good 
agreement between DFT-MD results and the SC 
model is found up to 50 000 K. For lower densi- 
ties, the SC model predicts a sharp entropy in- 
crease at 20 000 K, which is again a result of the 
treatment of ionization effects. This trend is not 
confirmed by the DFT-MD simulations. One also 
finds significant deviations at lower temperatures, 
in particular around rg=l.6. Even at a relatively 
low density of = 2.2, the agreement is not per- 
fect. From 4000 to 20 000 K, the SC model un- 
derestimates the entropy and it overestimates the 
entropy for lower temperatures. Figure [9] shows 
that such deviations persist over a wider density 
range. In principle, one expects a non- or weakly 




2000 5000 10000 20000 50000 10" 



Temperature (K) 

Fig. 6. — Helmholtz free energy per electron at 
constant density computed with DFT-MD simula- 
tions are compared with the analytical SC model. 
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Fig. 7. — Helmholtz free energy per electron at 
constant temperature computed with DFT-MD 
simulations are compared with the analytical SC 
model. 
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Fig. 8. — Entropy per electron at constant den- 
sity computed with DFT-MD simulations are com- 
pared with the analytical SC model. 
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interacting gas of hydrogen molecules and helium 
atoms to be perfectly described by the SC model. 
However, the density that we can efficiently study 
with DFT-MD simulations does not yet appear to 
be low enough for the deviations to decay to zero. 



DFT-MD 
SC model 
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Fig. 9. — Entropy per electron at constant tem- 
perature computed with DFT-MD simulations are 
compared with the analytical SC model. The up- 
per panel shows results over a wide temperature- 
density range while the lower panel zooms in on 
the low density regime where hydrogen occurs in 
molecular form when the temperature is below 
5000 K. 



4. Free Energy Fit for the Equation of 
State 

We fitted our ab initio results for P, E, F, 
and S in table [T] with a two-dimensional spline 
function that represents the Helmholtz free energy 
in terms of temperature, T, and electron density, 
n — Ng/V. By construction, this fit is thermody- 
namically consistent. We employ the same func- 
tional form that we used to represent the free en- 
ergy of hot, dense helium in (iMilitzerii200S) . ex- 
cept the splines here are functions of n rather than 
log(n). Table provides the free energy as well as 
the required derivatives on a number of (n, T) knot 
points. Atomic units are used throughout. 

To evaluate the fit for (n*,T*), we first con- 
struct a separate one-dimensional cubic spline 
function, Fn{T), for every density on a grid rang- 
ing from r5=3.581 to 0.536 (0.0670-20.0gcm-3). 
At every density, the free energy is given on a num- 
ber of temperature knots and its first derivative, 
^ I ^ , is specified for the highest and lowest tem- 
peratures. We construct a similar one-dimensional 
spline function that represents ^[^(T) at the 
smallest and largest density. We then evaluate 
all these splines functions at T* and construct a 
one-dimensional spline function Ft-' (n) from the 
free energy values and its first derivatives at the 
boundaries. This provides us not only with a 
straightforward way to obtain the free energy at 
every (n, T) point but we can also derive the pres- 
sure and entropy by taking analytical derivatives, 



P 



dF 
dn 



and S 



dF 
df 



(3) 



The internal energy and Gibbs free energy then 
follow from E ^ F + TS and G ^ F + PV. 

When we constructed this fit, we made sure 
every EOS point in table [T] is well reproduced. 
We extended the domain of the fit a bit beyond 
the range of the DFT-MD data. This leads to a 
smoother representation of the data in the interior 
of the domain and also allows us to gradually ap- 
proach the SC EOS in the limit of low density. As 
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Fig. 10. — Adiabats derived from DFT-MD simu- 
lations are compared with the SC model. The la- 
bels denote the entropy in units of kb per electron. 
The circles indicate parameters where entropy and 
free energy have been calculated in addition to the 
pressure and internal energy. The horizontal ar- 
rows label conditions where the deviation from the 
SC model are large and important for the interiors 
of Saturn and Jupiter. The vertical arrows indi- 
cate deviations are high temperature where the SC 
model does not treat electronic excitations accu- 
rately. 



figure [10] shows, we were able to smoothly match 
onto the SC adiabats for entropy values from 6 to 
10 and again for 13 and 14 fcb/el. A disagreement 
remains for S=ll and 12 fcb/el. but the SC EOS 
is not thermodynamically consistent in the regime 
of 10 000 to 20 000 K and no attempt was made 
to reproduce those adia bats. So far, only the one 
exoplanet HAT-P-32b (jHartman et al.ll201ll ) ap- 
pears to have an internal entropy in excess of 11 
fcfc/el.; see Figure [T^ 

We also find a significant disagreement in the 
high density limit between our DFT-MD adiabats 
and the predictions of the SC model. Starting with 
entropy values of S" = 10fcf,/el., the DFT-MD re- 
sults predict the adiabats reach states of higher 
temperatures and higher pressure for a given den- 
sity. 

The most significant result of Fig. [10] is the de- 
viations along the S—1 ki,/e\. adiabat. The DFT- 
MD simulations predict a decrease in slope of the 
adiabat exactly where the hy drogen molecules dis- 



sociate (jMilitzer et al.ll2008[ ). Since the SC model 



interpolates between separate atomic/metallic and 
molecular thermodynamic descriptions, it has no 
predictive power in the regime of pressure dis- 
sociation where all the different species interact 
strongly. 

5. Giant Planet Interiors 

In Fig. 1111 we compare different predictions 
for the adiabat i n Ju piter's interior. Similar 
to Militzer et al. ( 20081) . the calculations of the 
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Fig. 11. — Adiabats from different calculations for 
Jupiter's interior. 
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entropy by iNettelmann et aL (I2OO8L I2OI2I) relied 
solely on the P and E from DFT-MD simulations. 
Since no TDI was employed, the entropy was de- 
termined indirectly from thermodynamic relation- 
ships and an integration over a large path through 
density-temperature space. Here one faces two 
challenges. Since the integration can only yield 
the entropy difference between two p-T points, one 
needs to find a starting point for the integration 
where DFT-MD simulations work and the entropy 
is known reliably through other means. Secondly, 
one needs to determine P and £^ on a very fine 
grid in p-T space, so that integration errors do not 
accumulate. Since both challenges are difficult to 
meet, we instead adopted the more reliable TDI 
te chnique that w e extended to molecular systems 



Militzeil (|2013f ) 



Figure [TT] shows that there exist some dis- 
cre pancies between our TDI ca lculations and 
the INettelmann et all (|2008l I2OI2I ) results in the 
low density as well as the high density limit. In the 
molecular regime at 10 GPa, Nettelmann et al.'s 
temperatures for Jupiter's adiabat are 8% higher 
than our TDI calculations predict. In the metallic 
regime at 4000 GPa, their results are 19% higher 
than our TDI predictions. This implies that there 
exist discrepancies in the low and high density 
limits before any ab initio P and E data points 
ar e entered into the calculation of the adiabats 
(2008, 
)n the 



bv INettelmann et al 
Adiabats based 



INettelmann et al 



(j2012[ ) work now show a pronounced flattening 
in the regime of molecular disso ciation from 15-40 



GPa t hat was not present in the INettelmann et al 
( 20081 ) calculations. The pressure range simi- 
lar to our ab initio results but the magnitude 
is higher than we predict based on TDI. F rom 40 



to 200 GPa, the INettelmann et all (|2012D calcu 



lations predict a steep rise in temperature for 
the adiabats. This is not consistent with our 
TDI calculations and im plies that in models by 
Nettelmann et al.l (|2012| ) most of Jupiter's mass is 
at 19% higher temperature than we predict based 
on our TDI calculations. In the low and high den- 
sity limits, our adiabats are instead in relatively 
good agreement with the SC model. 



6. Mass-Radius Relationships 

The vexing problem of radius anoma l ies o f 
transiting giant planets ( Burrows et al.l l2007l ) 
has continued with t he addition of more objects 
(jLaughlin et al.l 1201 ll ). Figure [T2l shows measure- 
men ts posted in the online Exoplanet Encyclopae- 
dia ([Schneider et al.l 120111 ) as of late 2012 and 
different theoretical curves that we will discuss 
below. Briefly, the problem arises from a sig- 
niflcant population of exoplanets that have radii 
too large to be explained by thermal distention 
from retained primordial heat, and there is a fur- 
ther population with radii well below those ex- 
pected for primarily H-He composition even at 
zero temperature. Any point that falls below 
the S = G/ch/el. curve can be explained by in- 
voking the presence of a rocky core and/or ad- 
mixture of heavy ele ments, which reduces t he ra- 
dius for given mass dMiller fc Fortnevll201ll ). But 
it is not so simple to classify the population of 
anomalously distended giant exoplanets, for the 
degree of distention depends on such factors as 
the planet's a ge and degree of irradiation from 
the host star ( Fortnev fc Nettelmann 2010l ). and 
possible additiona l heating mechanism s such as 
ohmic dissipation (|Batvgin et al.ll201ir) . 

In order to most clearly exhibit the differences 
in predicted radius between the DFT-MD simula- 
tions and analytical SG model, we model a planet 
of mass M as a H-He object of constant entropy S 
and fixed helium fraction oi Y — 0.245 with nei- 
ther a rocky core nor heavy element component in 
the gas envelope. Since we do not have DFT-MD 
simulation data at very low densities, we switch 
back to the SC model below 0.0670 gcm"-^, the 
lowest density of the free energy fit to our DFT- 
MD data. 

As is well known ( Ghandrasekhaij Il957 ) , for- 
mally such an object has a precisely defined ra- 
dius where the temperature T, mass density p, and 
pressure P simultaneously go to zero. In a real 
object, the ideal-gas outer layers cannot persist 
in such an isentropic state and instead the tem- 
perature reaches a finite limit set by the effective 
temperature for the radiation balance in the outer 
layers. The radius as measured by transit obser- 
vations also depends on sources of slant opacity in 
these outer layers. As explored bv iBurrows et al 



(|2007l ). the value of such a radius can vary by sev- 
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Fig. 12.— Radius R (in units of Rj 70000 
km vs. mass M (in units of Mj = Jupiter's 
mass). The solid data points are measurements 
of transiting exoplanets. The curves show the 
R{M) relation predicted from two EOS calcula- 
tions (solid curves are for DFT-MD simulations, 
dashed curves are for the SC model) for fixed en- 
tropies of 5 = 6, 6.8 (Saturn; heavy curve), 7 
(Jupiter), 8, 9, 10, 11, and 12 fcb/el. The three 
open data points denote Saturn, HD 209458b, and 
Jupiter from left to right. 



eral 0.1 R,j, depending on the atmospheric model. 
Adjusting the parameters controlling the atmo- 
sphere can move a model closer to agreement with 
objects of inflated radii, but sometimes a mis- 
match remains. A major result of the present 
paper is that differences in the EOS alone can 
also lead to radius changes of several 0.1 Be- 
cause we consider only strict-adiabatic models, the 
deviations that we point out are entirely due to 
differences in the EOS at high pressure. In Fig- 
ure I13i we exhibit these differences for various en- 
tropies. For entropy values up to 9 fcfc/el., the 
DFT-MD calculations consistently predict smaller 
planet radii than the SC model, which is a direct 
consequence of the density enhancement on the 
adiabats around 100 GPa illustrated in Figs. [10] 
and [11] Figure [10] also shows that the DFT-MD 
and SC adiabats for 10, 11, and 12 fcb/el. cross 
over in the density range from 0.15 to 0.7gcm~^. 
This is the reason why the DFT-MD calculation 
predict larger planet radii than the SC model for 
massive planets with M > 0.5Mj but significantly 
smaller radii for light planets. The deviations be- 
tween the DFT-MD and SC predictions in Fig. [13] 
reach values up to approximately 0.4 Jupiter radii. 

The exoplanet HD 209458b (middle open data 
point in Figure 1121) fortuitously falls near an en- 
tropy S w 9.5 fcb/el. and mass M w 0.7Afj where 
AR « 0. Nevertheless, the interior T-p and P-p 
profiles in figures [14] and [15] differ significantly be- 
tween the DFT-MD and SC EOSs. These figures 
also compare interior profiles for simplified (pure 
H-He mixtures on an adiabat) models of Jupiter 
and Saturn. 

Figure [16] shows differences in evolutionary be- 
havior of our simplified planetary models. This 
figure plots the value of the central temperature, 
^'central, VS. the Central density, Pcontrai, for a range 
of adiabats and masses. During the evolution of 
a planet of constant mass, its central density in- 
creases monotonically while its central tempera- 
ture exhibits a maximum. During the initial con- 
traction, the temperature in the center increases at 
first as the material is subjected to increasing pres- 
sure. When a degenerate interior state is reached, 
the contraction ceases and the whole planet starts 
to cool. According to DFT-MD simulations the 
maximum temperature reached is up to 10 000 K 
lower than predicted by the SC model. This devi- 
ation may have consequences for the evolution of 
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cores in giant planets that remain to be explored. 
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Fig. 13. — The ordinate is the radius difference 
Ai? = i?DFT-MD — -Rsc for a given adiabat. The 
heavy shaded curve is for a Saturn adiabat with 
5 = 6.8. At low masses (below about 0.1 Mj for 
low entropies) the difference goes to zero because 
there are no DFT-MD data at low densities. 
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7. Conclusions 

This paper provides an equation state table for 
hydrogen-helium mixtures in giant planet interiors 
that was derived from ab initio computer simula- 
tions. The combination with an efficient thermo- 
dynamic integration technique enabled us to cal- 
culate the entropy and free energy directly, in ad- 
dition to pressure and internal energy that follow 
from standard simulations. 

Our complete EOS table with 391 density- 
temperature points as well as a thermodynamially 
consistent free energy fit is included in this publi- 
cation so that our EOS can be easily incorporated 
in future models for giant planet interiors. 

We have identified significant deviations for the 
Saumon and Chabrier EOS models. The new 
DFT-MD EOS causes low-entropy giant-planet 
models {S < 8fcf,/el.) to shrink in comparison 
to SC models by up to 0.08 Jupiter radii. But 
for hot giant planets with mass exceeding 0.5 Afj 
and with interior entropy values in the range from 
10-12 fcf,/cl., the DFT-MD simulations predict 
significantly larger radii. The correction to the 
SC model reaches 0.4 Jupiter radii for the hottest 
planets. Thus, the revision suggests that some 
of the most infiated giant exoplanets are at lower 
entropies than was previously inferred. Our re- 
vision could ameliorate the "inflated giant exo- 
planet" discrepancy to some extent but perhaps 
not for HD209458b. The matter is to be revisited 
with detailed evolutionary calculations based on 
our revised EOS. 

This work has been supported by NASA and 
NSF. Computational resources at NCCS were 
used. 



Fig. 14. — Temperature-density profile for the in- 
terior of three specific objects. Solid curves are 
for DFT-MD and dashed curves are for SC. The 
enhanced density in the DFT-MD model (see ar- 
row) at pressures near 100 GPa produces a slightly 
smaller radius for both Jupiter and Saturn. Below 
a density of 0.0670 gcm""^, the SC model is used 
for both calculations. 
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Fig. 15. — Pressure-density profile for the plane- 
tary interiors shown in Fig. 1141 




Fig. 16. — The values of T and p at the center of 
a planets according to the DFT-MD (heavy sohd 
lines) and SC (dashed Unes) EOSs for five differ- 
ent planet masses corresponding to 0.3 Mj (Sat- 
urn, lowest curves), 0.7 Mj (HD 209458b), 1.0 Mj 
(Jupiter), 1.5 Mj, and 2.0 Mj (top curves). The 
thin solid lines refer to various DFT-MD adiabats 
with entropy values from 6 (lowest curve) to 12 
(steepest curve) in steps of 0.5 A:f,/el. The heavy 
solid line is for Saturn's entropy of 6.8. 
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Table 1 

Equation of state derived from DFT-MD simulations. 





Density 


Temperature 


Pressure 


Internal Energy 


Helinholtz Free 


Entropy 


(ao) 


(gcm^-^) 


(K) 


(GPa) 


(Ha/el) 


Energy (Ha/el) 


(fcb/el) 


0.70 


8.9658 


5000 


17713.9(3) 


0.69251(3) 


0.641237(7) 


3.238(3) 


0.80 


6.0064 


5000 


8170.3(4) 


0.41280(4) 


0.351520(11) 


3.870(3) 


0.90 


4.2185 


5000 


4064.5(3) 


0.24343(4) 


0.173088(24) 


4.443(4) 


1.00 


3.0753 


5000 


2141.5(2) 


0.13726(3) 


0.058946(16) 


4.946(3) 


1.10 


2.3105 


5000 


1180.6(2) 


0.06926(6) 


-0.016209(11) 


5.398(4) 


1.20 


1.7797 


5000 


675.0(1) 


0.02513(5) 


-0.066864(21) 


5.810(4) 


1.30 


1.3998 


5000 


398.4(1) 


-0.00370(5) 


-0.101600(16) 


6.183(4) 


1.40 


1.1207 


5000 


242.1(1) 


-0.02265(8) 


-0.12587(3) 


6.519(7) 


1.50 


0.9112 


5000 


151.8(2) 


-0.03527(6) 


-0.14310(4) 


6.810(6) 


1.60 


0.7508 


5000 


98.0(1) 


-0.04400(6) 


-0.15565(2) 


7.051(5) 


1.86 


0.4779 


5000 


38.2(1) 


-0.05755(7) 


-0.17565(3) 


7.459(7) 


2.00 


0.3844 


5000 


25.74(8) 


-0.06295(15) 


-0.18271(4) 


7.563(12) 


2.10 


0.3321 


5000 


19.97(9) 


-0.06627(13) 


-0.18655(4) 


7.596(11) 


2.20 


0.2888 


5000 


15.51(9) 


-0.0685(2) 


-0.19015(4) 


7.683(16) 


2.30 


0.2528 


5000 


12.24(7) 


-0.0702(2) 


-0.19312(5) 


7.765(14) 


2.40 


0.2225 


5000 


9.64(6) 


-0.0714(2) 


-0.19580(7) 


7.855(17) 



Note. — Table [T] is published in its entirety in the electronic edition of the Astrophysical Journal. 
A portion is shown here for guidance regarding its form and content. 



Table 2 

Coefficients of Free Energy Fit for the equation of state. 







T 


coefficient 




(a.u.) 


(a.u.) 


(a.u.) 


BF 

Ut 


1.75 


0.001013 


-8.67942x10"^ 


F 


1.75 


0.001013 


-9.23942x10"^ 


F 


1.75 


0.003131 


-9.76575x10"^ 


F 


1.75 


0.006512 


-1.12223x10"^ 


F 


1.75 


0.011911 


-1.44120x10"^ 


F 


1.75 


0.020533 


-2.07646x10"^ 


F 


1.75 


0.034302 


-3.23222x10"^ 


F 


1.75 


0.056290 


-5.29287x10"^ 


F 


1.75 


0.091403 


-8.92177x10"^ 


F 


1.75 


0.147476 


-1.53056 


F 


1.75 


0.237021 


-2.64345 


F 


1.75 


0.380018 


-4.58672 


OF 
OT 


1.75 


0.380018 


-1.41575x10+^ 



Note. — Table [2] is published in its entirety 
in the electronic edition of the Astrophysical 
Journal. A portion is shown here for guidance 
regarding its form and content. While this 
manuscript is still under review, a copy 
of the EOS fit as well as the computer in- 
terpolation program may be requested via 
email: militzer@berkeley.edu 
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